maindir<-getwd()

eventmeans<-TRUE
source("./psid_eventpreamble.R")
#empyet ==0 means you were ever on SSI or didn't work before shock.
otherdata<-rawdata[empyet==0 & educ<=2,]
rawdata<-rawdata[empyet==1 & educ<=2,]
rawdata[,wly:=wly/1000]
rawdata[,marnow:=married]

gc()
source("./sipp_eventpreamble.R")
setwd(maindir)

sipp[,newid:=interaction(ssuid,epppnum,spanel)]
sipp[ejbhrs1==-8,ejbhrs1:=20]
sipp[ejbhrs2==-8,ejbhrs2:=20]
sipp[ejbhrs1==-1,ejbhrs1:=NA]
sipp[ejbhrs2==-1,ejbhrs2:=NA]
sipp[is.na(ehrsall),ehrsall:=ifelse(!is.na(ejbhrs1),ejbhrs1,0)+ifelse(!is.na(ejbhrs2),ejbhrs2,0)]
sipp[,emp:=ehrsall>=30]
sipp[,c_agebin:=floor(tage/10)*10]
sipp[,post := time_limit >= 0]


sipp[,marnow:=as.numeric(married)]

sipp <- sipp[!is.na(marnow)&
             as.numeric(eeducate) <= 39 & 
             as.numeric(eeducate)!=-1,]

sipp[,edisabl:=edisabl=="Yes"]
sipp[,edisprev:=edisprev=="Yes"]
sipp[eeducate!=-1,hsgrad:=eeducate>=39]
sipp[,white:=erace=="White alone"]
sipp[,n:=1]
sipp[!is.na(edisabl)&!is.na(edisprev),health:=1]
sipp[edisabl==1 & !is.na(edisprev),health:=2]
sipp[edisprev==1,health:=3]
sipp[,ageold:=tage>=45]

sipp[,sp_tpearn_real:=sp_tpearn*12/1000]
sipp[is.na(sp_tpearn_real),sp_tpearn_real:=0]
sipp[,w_no_bus:=(thhtnw-pmax(thhbeq,0))/(1000*price)]
sipp[,thhtnw:=thhtnw/1000]


##HRS##
#already excluded women, highly educated, ever self-employed, ever on SSI alone, or who never work while healthy
event<-"workprob"
gvar<-"haspartner"
source("./hrs_eventpreamble.R")
setwd(maindir)
hrs<-hrsdata

hrs[,w_no_bus:=(wealth_tot-wealth_bus)/1000]
hrs[,w_no_bus_real:=(wealth_tot_real-wealth_bus_real)/1000]
hrs[,marnow:=as.numeric(prepartner)]
hrs[,disab:=as.numeric(disab)-1]
hrs[,health:=disab+1]
hrs[,newid:=id]
hrs[,n:=1]
for(v in names(hrs)[grepl("health_",names(hrs)) | grepl("functional_",names(hrs))]){
    hrs[,eval(v):=get(v)%in%c("1.yes","2.can't do",1:3)]
}
hrs[,diagnostic:=health_cancr==1 | health_heart == 1 | health_lung==1 | health_strok == 1]
hrs[,functional:=functional_walkr==1 | functional_dress==1 | functional_bath==1 | functional_eat==1 | functional_bed==1 | functional_toilt==1]
hrs[,functional_count:= (functional_walkr==1) + (functional_dress==1) + (functional_bath==1) + (functional_eat==1) + (functional_bed==1) + (functional_toilt==1)]
for(v in c("health_psych","health_backprob","diagnostic","functional","functional_count")){
hrs[,eval(v):=as.numeric(get(v))]  
}
hrs[is.na(sp_ly_real),sp_ly_real:=0]
hrs[,sp_ly_real:=sp_ly_real/1000]


decimals<-2

stat_row<-function(rawdata,varname,variable,dec,weights=1,se=FALSE){
  coeftable<-feols(as.formula(paste0(variable, "~ as.factor(marnow) - 1")), data = rawdata, weights=weights,cluster="newid")$coeftable
  diff<-feols(as.formula(paste0(variable, "~ marnow")), data = rawdata, weights=weights,cluster="newid")$coeftable
  if(se==FALSE) {
    object <- TR(varname)%:%TR(coeftable[,"Estimate"],dec=dec#,
                   #pvalues=ifelse(is.na(coeftable[,"Pr(>|t|)"]),1,coeftable[,"Pr(>|t|)"])
                   )#%:%
#    TR(diff["marnow","Estimate"],dec=dec,
#       pvalues=ifelse(is.na(diff["marnow","Pr(>|t|)"]),1,diff["marnow","Pr(>|t|)"]))
  } 
  if(se == TRUE){
    object <- TR("")%:%TR(coeftable[,"Std. Error"],se=TRUE,dec=dec)#%:%
#    TR(diff["marnow","Std. Error"],se=TRUE,dec=dec)
  }
  return(object)
}

stat_addcol<-function(rawdata,variable,dec,weights=1,se=FALSE){
  coeftable<-feols(as.formula(paste0(variable, "~ as.factor(marnow) - 1")), data = rawdata, weights=weights,cluster="newid")$coeftable
  diff<-feols(as.formula(paste0(variable, "~ marnow")), data = rawdata, weights=weights,cluster="newid")$coeftable
  if(se==FALSE){
  object <- TR(c(coeftable[,"Estimate"]#,diff["marnow","Estimate"]
                 ),dec=dec#,
                   #pvalues=c(coeftable[,"Pr(>|t|)"]
                             #,diff["marnow","Pr(>|t|)"]
                  #           )
               )
  }
  if(se == TRUE){
    object <- TR(c(coeftable[,"Std. Error"]
                   #,diff["marnow","Std. Error"]
                   ),se=TRUE,dec=dec)
    }
  return(object)
}


tab<-
  TR(c("","(1)","(2)","(3)","(4)","(5)","(6)")) +
  TR(c("","PSID","SIPP","HRS"),cspan=c(1,2,2,2))+
  TR(c("","Single","Partnered","Single","Partnered","Single","Partnered"))+
  midrulep(list(c(2,3),c(4,5),c(6,7)))  +
  TR("Panel A: Individual Characteristics",cspan=4, surround = "\\textit{%s}")+
  stat_row(rawdata,"Age at Onset","onset_age",2)%:%stat_addcol(sipp,"age_limit",2)%:%stat_addcol(hrs,"refonset_age",2)+
  stat_row(rawdata,"Age at Onset","onset_age",2,se=TRUE)%:%stat_addcol(sipp,"age_limit",2,se=TRUE)%:%stat_addcol(hrs,"refonset_age",2,se=TRUE)+
  stat_row(rawdata,"High School Grad","hsgrad",2)%:%stat_addcol(sipp,"hsgrad",2)%:%stat_addcol(hrs,"hsgrad",2)+
  stat_row(rawdata,"High School Grad","hsgrad",2,se=TRUE)%:%stat_addcol(sipp,"hsgrad",2,se=TRUE)%:%stat_addcol(hrs,"hsgrad",2,se=TRUE)+
  stat_row(rawdata,"White","white",2)%:%stat_addcol(sipp,"white",2)%:%stat_addcol(hrs,"white",2)+
  stat_row(rawdata,"White","white",2,se=TRUE)%:%stat_addcol(sipp,"white",2,se=TRUE)%:%stat_addcol(hrs,"white",2,se=TRUE)+
  stat_row(rawdata,"Spousal Annual Income","wly",2)%:%stat_addcol(sipp,"sp_tpearn_real",2)%:%stat_addcol(hrs,"sp_ly_real",2)+
  stat_row(rawdata,"Spousal Annual Income","wly",2,se=TRUE)%:%stat_addcol(sipp,"sp_tpearn_real",2,se=TRUE)%:%stat_addcol(hrs,"sp_ly_real",2,se=TRUE)+
  stat_row(rawdata[time<0,],"Pre-Disability Net Wealth","w_no_bus",2)%:%stat_addcol(sipp[time_limit<0,],"w_no_bus",2)%:%stat_addcol(hrs,"w_no_bus_real",2)+
  stat_row(rawdata[time<0,],"Pre-Disability Net Wealth","w_no_bus",2,se=TRUE)%:%stat_addcol(sipp[time_limit<0,],"w_no_bus",2,se=TRUE)%:%stat_addcol(hrs,"w_no_bus_real",2,se=TRUE)+
  vspace(5)+ 
  TR("Panel B: Panel Characteristics",cspan=4, surround="\\textit{%s}")+
  TR("Total Households")%:%TR(as.character(c(unique(rawdata[,.(newid,marnow)])[,.N,by=marnow][order(marnow),]$N)))%:% 
                              TR(as.character(c(unique(sipp[,.(newid,marnow)])[,.N,by=marnow][order(marnow),]$N)))%:% 
                              TR(as.character(c(unique(hrs[,.(newid,marnow)])[,.N,by=marnow][order(marnow),]$N)))+ 
  TR("Years Observed per Household")%:%TR(c(rawdata[,sum(n),by=.(newid,marnow)][,mean(V1),by=marnow][order(marnow),]$V1),dec=1)%:% 
                              TR(c(unique(sipp[,.(newid,rhcalyr,marnow,n)])[,sum(n),by=.(newid,marnow)][,mean(V1),by=.(marnow)][order(marnow),]$V1),dec=1)%:% 
                              TR(c(hrs[,sum(n),by=.(newid,marnow)][,mean(V1),by=marnow][order(marnow),]$V1),dec=1) 


TS(tab,file="all_stats",
   output_path=".",
   pretty_rules=T,
   header=c('r',rep('c',6)))

  
  
  
  ###########################################
  #Diagnostic measures among severely disabled only:#
  ###########################################

  v<-"marnow"
  
  vref <- 1
  outnum<-1
  tabhead<-  TR(c("","(1)","(2)","(3)")) +
    TR(c("","Single","Partnered","Partnered - Single"),cspan=c(1,1,1))+
    midrulep(list(c(2,4)))  
  #    vspace(5)
  taba<- TR(c("Panel A: PSID Sample"),cspan=c(4), surround = "\\textit{%s}")
  
  outnames<-c("Suffered Cancer/Cardiovascular Event","Any Functional Limitation",
              "Had Blue Collar Job when Last Healthy")
  #"Total Functional Limitations")
  for(out in c("diagnostic","trouble","bluecollar")){
    #,"trouble_count"
    rawdata[,pweight:=NULL]
    rawdata[,pval:=NULL]
    rawdata[!is.na(get(v)) & !is.na(get(out)) & !is.na(agebin) &health==3
            ,pval:=feols(as.formula(paste0(v,"~ 1 | interaction(as.factor(agebin),as.factor(health))")),
                         data = rawdata[!is.na(get(v)) & !is.na(get(out)) & !is.na(agebin)  &health==3,
                         ])$fitted.values]
    
    rawdata[get(v)==vref & !is.na(pval),pweight:=1]
    rawdata[get(v)==0,
            pweight:=pval/(1-pval)]
    #Imposing common support:
    rawdata[(round(pval,6) == 1 | round(pval,6) == 0),pweight:=NA]
    rawdata[!is.na(pweight)&!is.na(get(v)),weighted.mean(age,pweight),by=c(eval(v),"health")]
    rawdata[,mean(age),by=c(eval(v),"health")]
    #######################################
    
    means<-feols(as.formula(paste0(out,"~interaction(as.factor(",v,'),as.factor(health))-1')),
                 data= rawdata[,],
                 weights=rawdata[,pweight],
                 cluster="newid")
    diffs<-feols(as.formula(paste0(out,"~as.factor(health) + as.factor(health):",v)),
                 data= rawdata[,],
                 weights=rawdata[,pweight],
                 cluster="newid")
    taba<-taba+
      TR(c(outnames[outnum]))%:%TR(c(means$coefficients[grepl("0.",names(means$coefficients))], means$coefficients[grepl("1.",names(means$coefficients))]
      ),dec=decimals,
      pvalues=c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"], means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"]
      ),
      )%:%TR(c(diffs$coefficients[grepl("marnow",names(diffs$coefficients))]),dec=decimals,
             pvalues=diffs$coeftable[grepl("marnow",names(diffs$coefficients)),"Pr(>|t|)"])+
      TR(c(""))%:%TR(c(means$se[grepl("0.",names(means$se))], means$se[grepl("1.",names(means$se))]
      ),dec=decimals,se=TRUE)%:%TR(c(diffs$se[grepl("marnow",names(diffs$se))]),dec=decimals,se=TRUE)
    outnum<-outnum+1
  }
  

  outnum<-1
  tabc<-TR(c("Panel B: SIPP Sample"),cspan=4, surround = "\\textit{%s}")
  outnames<-c("Any Functional Limitation","Total Functional Limitations","Limitation: Mental Disorder","Limitation: Congenital Disorder","Limitation: Musculoskeletal Disorder","Limitation: Cancer/Cardiovascular Event")
  for(out in c("functional","functional_count","mental","congenital","musculoskeletal","diagnostic")){
    sipp[,pweight:=NULL]
    sipp[,pval:=NULL]
    sipp[!is.na(get(v)) & !is.na(get(out)) & !is.na(c_agebin) &health==3
         ,pval:=feols(as.formula(paste0(v,"~ 1 | interaction(as.factor(c_agebin),as.factor(health))")),
                      data = sipp[!is.na(get(v)) & !is.na(get(out)) & !is.na(c_agebin)  &health==3,
                      ])$fitted.values]
    
    sipp[get(v)==vref & !is.na(pval),pweight:=1]
    sipp[get(v)==0,
         pweight:=pval/(1-pval)]
    #Imposing common support:
    sipp[(round(pval,6) == 1 | round(pval,6) == 0),pweight:=NA]
    sipp[!is.na(pweight)&!is.na(get(v)),weighted.mean(tage,pweight),by=c(eval(v),"health")]
    sipp[,mean(tage),by=c(eval(v),"health")]
    #######################################
    
    print(paste0("Unweighted, Condition: ",out))
    print(feols(as.formula(paste0(out,"~as.factor(",v,')')),
                data= sipp,
                cluster="ssuid")
    )
    
    print(paste0("Weighted, Condition: ",out))
    print(feols(as.formula(paste0(out,"~as.factor(",v,')')),
                data= sipp,
                weights=sipp[,pweight],
                cluster="ssuid")
    )
    means<-feols(as.formula(paste0(out,"~interaction(as.factor(",v,'),as.factor(health))-1')),
                 data= sipp,
                 weights=sipp[,pweight],
                 cluster="ssuid")
    diffs<-feols(as.formula(paste0(out,"~as.factor(health) + as.factor(health):",v)),
                 data= sipp[,],
                 weights=sipp[,pweight],
                 cluster="ssuid")
    tabc<-tabc+
      TR(c(outnames[outnum]))%:%TR(c(means$coefficients[grepl("0.",names(means$coefficients))],means$coefficients[grepl("1.",names(means$coefficients))]
      ),dec=decimals,
      pvalues=ifelse(is.na(c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"], means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"]
      )),
      1,
      c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"], means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"]
      )),
      )%:%TR(c(diffs$coefficients[grepl("marnow",names(diffs$coefficients))]),dec=decimals,
             pvalues=ifelse(is.na(diffs$coeftable[grepl("marnow",names(diffs$coefficients)),"Pr(>|t|)"]),
                            1,
                            diffs$coeftable[grepl("marnow",names(diffs$coefficients)),"Pr(>|t|)"]))+
      TR(c(""))%:%TR(c(means$se[grepl("0.",names(means$se))] , means$se[grepl("1.",names(means$se))]
      ),dec=decimals,se=TRUE)%:%TR(c(diffs$se[grepl("marnow",names(diffs$se))]),dec=decimals,se=TRUE)
    outnum<-outnum+1
  }
  
  outnum<-1
  tabc_workchecked<-TR(c("Panel C: SIPP Sample, Observed Working"),cspan=4, surround = "\\textit{%s}")
  outnames<-c("Any Functional Limitation","Total Functional Limitations","Limitation: Mental Disorder","Limitation: Congenital Disorder","Limitation: Musculoskeletal Disorder","Limitation: Cancer/Cardiovascular Event")
  for(out in c("functional","functional_count","mental","congenital","musculoskeletal","diagnostic")){
    sipp[,pweight:=NULL]
    sipp[,pval:=NULL]
    sipp[!is.na(get(v)) & !is.na(get(out)) & workcheck==1 & !is.na(c_agebin) &health==3
         ,pval:=feols(as.formula(paste0(v,"~ 1 | interaction(as.factor(c_agebin),as.factor(health))")),
                      data = sipp[!is.na(get(v)) & !is.na(get(out)) & workcheck==1 & !is.na(c_agebin)  &health==3,
                      ])$fitted.values]
    
    sipp[get(v)==vref & !is.na(pval),pweight:=1]
    sipp[get(v)==0,
         pweight:=pval/(1-pval)]
    #Imposing common support:
    sipp[(round(pval,6) == 1 | round(pval,6) == 0),pweight:=NA]
    sipp[!is.na(pweight)&!is.na(get(v)),weighted.mean(tage,pweight),by=c(eval(v),"health")]
    sipp[,mean(tage),by=c(eval(v),"health")]
    #######################################
    
    print(paste0("Unweighted, Condition: ",out))
    print(feols(as.formula(paste0(out,"~as.factor(",v,')')),
                data= sipp,
                cluster="ssuid")
    )
    
    print(paste0("Weighted, Condition: ",out))
    print(feols(as.formula(paste0(out,"~as.factor(",v,')')),
                data= sipp,
                weights=sipp[,pweight],
                cluster="ssuid")
    )
    means<-feols(as.formula(paste0(out,"~interaction(as.factor(",v,'),as.factor(health))-1')),
                 data= sipp,
                 weights=sipp[,pweight],
                 cluster="ssuid")
    diffs<-feols(as.formula(paste0(out,"~as.factor(health) + as.factor(health):",v)),
                 data= sipp[,],
                 weights=sipp[,pweight],
                 cluster="ssuid")
    tabc_workchecked<-tabc_workchecked+
      TR(c(outnames[outnum]))%:%TR(c(means$coefficients[grepl("0.",names(means$coefficients))],means$coefficients[grepl("1.",names(means$coefficients))]
      ),dec=decimals,
      pvalues=ifelse(is.na(c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"], means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"]
      )),
      1,
      c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"], means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"]
      )),
      )%:%TR(c(diffs$coefficients[grepl("marnow",names(diffs$coefficients))]),dec=decimals,
             pvalues=ifelse(is.na(diffs$coeftable[grepl("marnow",names(diffs$coefficients)),"Pr(>|t|)"]),
                            1,
                            diffs$coeftable[grepl("marnow",names(diffs$coefficients)),"Pr(>|t|)"]))+
      TR(c(""))%:%TR(c(means$se[grepl("0.",names(means$se))] , means$se[grepl("1.",names(means$se))]
      ),dec=decimals,se=TRUE)%:%TR(c(diffs$se[grepl("marnow",names(diffs$se))]),dec=decimals,se=TRUE)
    outnum<-outnum+1
  }
  
  outnum<-1
  tabd<-TR(c("Panel C: HRS Sample"),cspan=4, surround = "\\textit{%s}")
  outnames<-c("Suffered Mental Disorder","Suffered Musculoskeletal Disorder","Suffered Cancer/Cardiovascular Event","Any Functional Limitation","Total Functional Limitations")
  for(out in c("health_psych","health_backprob","diagnostic","functional","functional_count")){
    hrs[,pweight:=NULL]
    hrs[,pval:=NULL]
    hrs[!is.na(get(v)) & !is.na(get(out)) & !is.na(agebin) & health == 2
        ,pval:=feols(as.formula(paste0(v,"~ 1 | interaction(as.factor(agebin),as.factor(health))")),
                     data = hrs[!is.na(get(v)) & !is.na(get(out)) & !is.na(agebin)  & health == 2,
                     ])$fitted.values]
    
    hrs[get(v)==vref & !is.na(pval),pweight:=1]
    hrs[get(v)==0,
        pweight:=pval/(1-pval)]
    #Imposing common support:
    hrs[(round(pval,6) == 1 | round(pval,4) == 0),pweight:=NA]
    hrs[!is.na(pweight)&!is.na(get(v)),weighted.mean(age,pweight),by=c(eval(v),"health")]
    hrs[,mean(age),by=c(eval(v),"health")]
    #######################################
    
    means<-feols(as.formula(paste0(out,"~interaction(as.factor(",v,'),as.factor(health))-1')),
                 data= hrs,
                 weights=hrs[,pweight],
                 cluster="newid")
    diffs<-feols(as.formula(paste0(out,"~as.factor(health) + as.factor(health):",v)),
                 data= hrs[,],
                 weights=hrs[,pweight],
                 cluster="newid")
    tabd<-tabd+
      TR(c(outnames[outnum]))%:%TR(c(means$coefficients[grepl("0.",names(means$coefficients))]),dec=decimals,
                                   pvalues=ifelse(is.na(c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"])),
                                                  1,
                                                  c(means$coeftable[grepl("0.",names(means$coefficients)),"Pr(>|t|)"])))%:%
      TR(c(means$coefficients[grepl("1.",names(means$coefficients))]),dec=decimals,
                                  pvalues=ifelse(is.na(c(means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"])),
                                                 1,
                                                 c(means$coeftable[grepl("1.",names(means$coefficients)),"Pr(>|t|)"])))%:%
      TR(c(diffs$coefficients[grepl("marnow",names(diffs$coefficients))]),dec=decimals,
         pvalues=diffs$coeftable[grepl("marnow",names(diffs$coefficients)),"Pr(>|t|)"])+
      TR(c(""))%:%
      TR(c(means$se[grepl("0.",names(means$se))]),dec=decimals,se=TRUE)%:%
      TR(c(means$se[grepl("1.",names(means$se))]),dec=decimals,se=TRUE)%:%
      TR(c(diffs$se[grepl("marnow",names(diffs$se))]),dec=decimals,se=TRUE)
    outnum<-outnum+1
  }
  

  TS(tabhead +
       tabc,file=paste0(v,"_sev_table_c"),
     output_path=".",
     pretty_rules=T,
     header=c('r',rep('c',3)))
  
  TS(tabhead +
       tabc_workchecked,file=paste0(v,"_sev_table_c_workchecked"),
     output_path=".",
     pretty_rules=T,
     header=c('r',rep('c',3)))
  
  TS(tabhead +
       tabd,file=paste0(v,"_sev_table_d"),
     output_path=".",
     pretty_rules=T,
     header=c('r',rep('c',3)))
  
  
  TS(tabhead + vspace(5)+
       taba+ vspace(10)+
       #       tabb+ vspace(10)+
       tabc+ vspace(10)+
       tabd,file=paste0(v,"_sev_table"),
     output_path=".",
     pretty_rules=T,
     header=c('r',rep('c',3)))
  
